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We study the dynamical behavior of disordered many-particle systems with long-range Coulomb 
interactions by means of damage-spreading simulations. In this type of Monte-Carlo simulations 
one investigates the time evolution of the damage, i.e. the difference of the occupation numbers of 
two systems, subjected to the same thermal noise. We analyze the dependence of the damage on 
temperature and disorder strength. For zero disorder the spreading transition coincides with the 
equilibrium phase transition, whereas for finite disorder, we find evidence for a dynamical phase 
transition well below the transition temperature of the pure system. 
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O ■ I. INTRODUCTION 

The combined influence of disorder and long-range interactions on the properties of many-particle systems has been 
a subject of great interest for some time. In electronic systems already disorder or interactions alone can drastically 
change the physical behavior. Disorder can lead, e.g., to a metal-insulator transition due to Anderson localization. On 
c*3 i the other hand, a metal-insulator transition can also be induced by correlations due to electron-electron interactions. 
If disorder and interactions are both significant then complex physical problems and phenomena arise, many of which 

S are not completely understood. 
The behavior of strongly localized correlated electrons in disordered insulators is especially complicated, both 
experimentally and theoretically. Thus progress has been slow since the first investigations.™] Many properties of 
such systems are still poorly understood. In particular there are only few and contradicting results on thermodynamics, 
phase diasram, phase transitions or critical behavior, and the examination of the dynamical behavior is only at its 
O , beginning.!] Two of the central questions are whether or not the disordered interacting electron system shows glassy 
behavior and what is the nature of the glassy "state". Two different views can be found in the literature. In the 
t—{ \ earlier work the formal similarity between disordered localized electrons and spin glasses had lead to speculations 
■ about a possible equilibrium phase transition to a spin-glass-like low-temperature phase.uu More recent investigations 
0^ ' show, however, growing experimental and theoretical evidence of the transition being of dynamical natureHu 

In this paper we study the dynamical behavior of disordered localized electrons by means of the damage-spreading 
method. In this type of Monte-Carlo simulations the microscopic differences of the time evolution between two systems 
are investigated. In particular, we address the question of a dynamical phase transitions from a dynamically active 
high-temperature phase to a frozen low-temperature phase upon changing characteristic parameters like disorder or 
temperature. Our paper is organized as follows. In Sect. || we introduce the Coulomb glass model, the prototype 
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model of disordered localized electrons. In Sect. Ill we describe the damage spreading technique, whereas in Sect. 
IV we present the results for the dynamical behavior of the model. Section [v| is dedicated to some discussions and 



conclusions. 



II. MODEL 



Our investigations are based on the Coulomb glass model first proposed by Efros and Shklavskiio to describe 
compensatfid-doped semiconductors. Later it has also been applied to simulate granular metalsE2l and conducting 
polymersOll3 The model consists of a square or cubic lattice of linear size L with N = L d sites (in d dimensions) and 
lattice constant a. The sites can be occupied by KN (0 < K < 1) electrons. These electrons are interacting via an 
unscreened Coulomb potential. To guarantee charge neutrality every site carries a compensating charge of +Ke (— e 
is the charge of the electron). The disorder of this system is described by the random potential The Hamiltonian 
of the Coulomb glass is given by 

H = Y,{y i -li)n i + \Y,{n i -K){n j -K)U ij = — (1) 

where fi is the chemical potential, rii (with values or 1) is the occupation number of site i and ry denotes the 
distance between sites i and j. In the rest of the paper we set the interaction strength between nearest neighbor 
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sites e 2 /a = 1 which fixes the energy scale. The random potential energies ipi are independent from each other and 
chosen according to some probability distribution W(cp). We use the box distribution with mean and width Wq. 
The parameter Wo measures the strength of the disorder. Specifically, we investigate a half-filled system (K = I/2). 
Then the Coulomb glass model is particle-hole symmetric and the chemical potential vanishes. (Note that the two 
quantities K and fi are not independent of each other. We treat K as a free parameter and calculate [i from it.) 

For later reference we briefly mention some properties of the Coulomb glass model. One of the central quantities is 
the single-electron density of states 



(2) 



at energy e and temperature T, where (. . .) denotes thermal and disorder averages. £j are the single-electron energies 
given by 



£i = (fii - (J, + ^ U ij( n j ~ K )- 



(3) 



The single-electron density of states of the Coulomb glass shows a pronounced gap, called the Coulomb gap, close 
to the Fermi energy ef (see Fig. l|). At zero temperature the density of states actually vanishes at the Fermi energyo, 
close to the Fermi energy it can be described by a power law 

g(e)<x\e-e F \ a (4) 

where a is approximately 1.2 for two-dimensional (2D) and 2.5 fop-3D systems. El At finite temperatures the Coulomb 
gap is filled gradually (for recent simulation results see, e.g., Ref.Ej). 
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FIG. 1. Single-electron density of states of the Coulomb glass at T = 0.008 for different strengths of disorder, em indicates 
the Madelung energy, ef the Fermi energy. 

The Coulomb glass model (Q) describes a system without internal dynamics. In reality the electrons, though 
localized, are coupled to additional (vibrational) degrees of freedom, which lead to transitions between the many- 
electron states. Phenomenologically this can be simulated by a Monte-Carlo method. In every Monte-Carlo step we 
change the occupation numbers of one or several sites with a certain probability. Within the Metropolis algorithm 
this probability is given by 



P = 



1 



,AH < 



exp[-M>] ,Aff>0 



(5) 
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where AH is the energy difference between the many-particle states before and after such a change, and ks is the 
Boltzmann constant. N such Monte-Carlo steps are called a Monte-Carlo sweep which is the natural time scale of 
our calculations. 

To simulate the dynamics one can use different "move classes" , which determine how the occupation numbers are 
changed in every Monte-Carlo step to get the new configuration. The simplest move class consists of exchanging 
a single electron with a reservoir (i.e. the conduction band in the case of doped semiconductors), other classes 
include hopping of single electrons between the sites, or correlated hopping of several electrons. In this paper we 
present results obtained by using only single-electron exchanges between the system and a reservoir, but we have also 
checked more complicated move classes. As long as we do not include distance-dependent "tunneling terms" into the 
transition probabilities (^), applying different move classes yields data which do not show a qualitatively different 
behavior. We attribute this result to the fact that single and multiple electron hops can be combined from the moves 
in our implementation of single electron exchanges with an external reservoir. Thus all many-electron states with 
K N electrons are available in our simulation. A more detailed investigation of this question including the effects of 
distance-dependent transition probabilities on the damage-spreading simulations is in progress. 



III. DAMAGE SPREADING 

The damage-spreading techniqueEl is a modification of the usual monte-Carlo method. The idea is to look not 
at the time evolution of a single system but to compare the time evolutions of two systems which are subjected to 
the same thermal noise (i.e., the same random numbers are used within the Metropolis algorithm). Usually, at the 
beginning of the simulation the occupation numbers of both systems differ only at a single site (or at a few sites, e.g. 
a single column in a 2D lattice system). 

Since both systems are thermodynamically identical, averages of equilibrium quantities will be the same for both 
systems. Microscopically, however, the two systems may evolve differently from each other. The central observable in 
damage-spreading simulations is the Hamming distance D(t), which is the portion of sites for which the occupation 
numbers differ between the two systems. D(t), which measures the "damage", is given by 

£»(*) = ^£K(t) ( fi ) 

i 

where n°(i)and n£(t) are the occupation numbers of site i of the original system and the copy at (Monte-Carlo) time 
t. For D(t) = the two systems are identical, D(t) = I/2 describes completely uncorrelated configurations, and for 
D{t) = 1 the two systems are totally anticorrelated. In the course of the time evolution the two systems evolve 
towards a steady state, in which D(t) fluctuates around an asymptotic average value 

1 f t+T 

D = lim lim - / dt'D{t') (7) 

r — >oo t — >oo T J j. 

Depending on the values of the external parameters temperature and disorder different regimes can be observed in 
principle if the initial damage D(0) is small: The damage may heal out during the time evolution (D = 0), the systems 
may stay partially correlated for infinite time (D < I/2), or the systems may become completely uncorrelated so that 
D = I/2. In contrast to the thermodynamics the detailed behavior of D(i) depends on the choice of the dynamical 
algorithm. Whereas Metropolis, Glauber and heat-bath dynamics give the same results for equilibrium quantities of 
a single system, the damage spreading results differ. For the Metropolis dynamics which we use (as well as foi^-thc 
Glauber dynamics) the damage tends to heal at low temperatures and tends to spread at high temperatures.E^I 
contrast, the heat-bath dynamics yields healing at high temperatures and frozen configurations at low temperatures. 
(Note that since D is not a thermodynamic quantity but measures the microscopic differences between two systems, 
there is no reason to expect that different dynamical algorithms give the same results.) 

We apply the damage-spreading technique to the 2D Coulomb-glass model at half filling K — 1/2 and linear system 
sizes L = 20, 80. The simulation proceeds as follows: (i) We create the initial system by choosing random potential 
values according to the probability distribution W(tp) and occupy the sites at random with KN electrons, (ii) We 
equilibrate this system at temperature T by performing several (at least 300) Monte-Carlo sweeps according to the 
Metropolis algorithm, (iii) A copy of the system is created and modified at a single site (or several sites). This 
difference in the occupation numbers constitutes the initial damage, (iv) We study the time evolution of the original 
and the copy using the same random numbers in the Metropolis algorithm for both systems. The damage D(t) is 
recorded and its asymptotic value D is determined. 
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Note that there is a modification of the, .damage-spreading method that can be used to determine equilibrium 
quantities instead of purely dynamic onesJij'Ej In that kind of simulations the occupation number of a single site in 
one of the systems is fixed whereas it is allowed to fluctuate in the other system. Consequently, the two systems are 
thcrmodynamically different and the damage can be related to equilibrium correlation functions. Since in this paper 
we are interested in the properties of the dynamics rather than in equilibrium quantities, our data is gained by means 
of the original damage-spreading method, where the occupation numbers of both systems are allowed to fluctuate. 

IV. RESULTS 
A. Time evolution 

In this subsection we present data on the time evolution of the dam age— D(t) starting with an initial damage 
consisting of a single site. In analogy to the well studied 2D Ising modelEjOEJ we find that for temperatures below 
a certain temperature T s , called the spreading temperature, the damage D(t) remains small and eventually heals, 
giving an asymptotic value of D — 0. For temperatures larger than T s the damage increases with time until a steady 
state is reached where D(t) fluctuates around a finite value. Consequently, the asymptotic damage D is finite in this 
regime. In Fig. || the time evolution of D(t) is shown for the Coulomb glass with zero disorder Wq = 0. 
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FIG. 2. Time dependence of the Hamming distance of the 2D Coulomb glass for different temperatures and Wo = 0. 



The three curves presented correspond to the three regimes discussed in the last section. At T = 0.5 the damage 
increases quickly and then fluctuates around D = lfe. This means the two systems become completely uncorrelated 
very fast. Consequently we are above the spreading temperature T s . At T = 0.1 the evolution of D(t) is much slower 
and the asymptotic damage is smaller than I/2. This behavior occurs, because the system is in the vicinity of the 
spreading transition at T s . It corresponds to the critical slowing down in ordinary critical phenomena. At T = 0.06 
the damage remains small and eventually heals, thus the system is below the spreading temperature T s . In the case 
of finite disorder Wq the time evolution of the damage is similar (see Fig. ||). 
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FIG. 3. Time dependence of the Hamming distance of the 2D Coulomb glass for different temperatures and Wo = 0.5. 



The asymptotic damage D is, however, different from or I/2 even far away from the spreading transition. The 
dependence of the damage on the external parameters temperature and disorder is investigated in more detail in 
Subsect. HVC. 



B. Influence of the long-range interaction 



The character of the interaction has a large influence on the time evolution of the damage. In systems with nearest- 
neighbor interactions, e.g. the Ising model, the damage can only spread within a single Monte-Carlo step from one 
site of the system to its neighbor. Therefore the clouds of damaged sites can grow only slowly in space and tend to be 
more compact (but not necessarily connected). In contrast, in systems with long-range interactions the occupation 
number of any site effects all other sites. The damage can spread from one site of the system to any other site within 
a single Monte-Carlo step. Therefore the damage spreads much faster as in systems with short-range interactions and 
the damage clouds are usually not compact. A comparison of the two cases is presented in Fig. |]. 




FIG. 4. Snapshot of the damage for 2D systems with short-range interactions (left) and long-range interactions (right) for 
T = 0.5 and Wo = 0.5 at a time of 5 Monte-Carlo sweeps after the introduction of a single damaged site. A filled circle indicates 
a damaged site where the occupation numbers of the two systems differ, an empty circle indicates that the occupation numbers 
of that site are identical in both systems. 
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Note, that since the damage can spread from one site to any other site in the case of long-range interactions, some of 
the methods developed to analyze the damage-spreading simulationstJO cannot be used for systems with long-range 
interactions. This applies to all methods that measure the spatial extension of the damage and its evolution, because 
the spatial extent of the damage cloud is not a well defined quantity for systems with long-range interactions. 

C. Temperature and disorder dependence of the asymptotic damage 

We now turn to the main results of this paper. Figure ^| shows an overview of the temperature and disorder 
dependence of the asymptotic Hamming distance D. 
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FIG. 5. Overview of the temperature dependence of the Hamming distance for various strengths of disorder of a 2D system 
with N = 20 2 sites. 

For disorder strength Wo = there is a pronounced transition at a spreading temperature of approximately T s = 0.1 
between a low-temperature regime with D = and a high-temperature regime with D = Ife. Within our numerical 
accuracy the spreading temperature T s coincides with the equilibrium critical point T c of the model without disorder 
which we determined from the peak in the specific heat C v of the Coulomb glass model as a function of temperature 
(see Fig. ||). 
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FIG. 6. Specific heat of the 2D Coulomb glass at Wo = calculated via the derivation of the internal energy of the system. 



For very high temperatures T — > 00 the spreading of the damage is drastically slowed down due to the fact that 
the probability P in the Metropolis algorithm, Eq. |5|, becomes independent of the actual configurations of the two 
systems (original and copy) and reaches P = 1 . This means that in both systems nearly every exchange of electrons is 
performed and differences in the occupation numbers occur only rarely. Our investigations of the spreading behavior 
for very high temperatures show that the Hamming distance D still reaches a plateau if plotted versus time as in 
Figs. and [| but the relaxation time diverges as is predicted in a recent mean-field theoryBiJ The damage-spreading 
transition in the Coulomb glass model without disorder occurs thus in complete analogy to that in the Ising model.Ej 

For finite disorder strength Wo, however, this behavior changes in several aspects. First, the values of the asymptotic 
Hamming distance in the high-temperature regime are smaller than D — I/2. This means, the two systems remain 
partially correlated even for high temperatures. The reason for that is easy to understand: In the presence of a random 
potential the electrons are trapped (repulsed) at sites with small (high) potential values ipi. These sites are identical 
in the original system and its copy. Therefore the presence of a random potential tends to reduce the damage. With 
increasing strength of disorder this trapping effect becomes larger, so that the maximum value of the damage is more 
and more reduced. On the other hand, increasing temperature makes it easier to overcome the potential differences 
so that the described reduction of the damage becomes less effective. 

The second effect of the disorder concerns the behavior of D at low temperatures and close to the spreading point. 
This region is shown in more detail in Fig. 
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In the case of finite disorder the asymptotic damage remains finite even at temperatures below the spreading 
temperature of the model without disorder. This somewhat counterintuitive result, viz. an acceleration of the 
dynamics by disorder, can be understood by looking at the single-electron density of states of the Coulomb glass 
model (see Fig. |l|). For Wo = the single-electron density of states at low temperatures has a hard gap around the 
Fermi energy Sp = and two peaks at the Madelung energies ±£m- Therefore there are only exponentially few sites 
that can be excited at low temperatures and thus the Hamming distance vanishes. In contrast, for finite disorder Wo, 
the gap in the density of states is not exponential but the power-law Coulomb gap ([|). Therefore more sites can be 
excited at low temperatures and the dynamics does not freeze completely, i.e., the Hamming distance remains finite. 

As can be seen in Fig. ^, even for finite disorder strength Wo there is, however, a spreading temperature T s (Wo), 
below which the asymptotic damage vanishes. T s (Wq) decreases with increasing Wo, but seems to tend to a finite 
limiting value for large Wq which we approximately determined to T s (oo) ~ 0.03. Note, that the existence of a 
spreading transition in the case of finite disorder is a purely dynamic phenomenon, since the system does not undergo 
an equilibrium phase transition. 
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FIG. 8. Hamming distance versus temperature for various system sizes at Wo = 0. 

In order to determine more detailed properties of the spreading transition a careful analysis of finite size effects is 
necessary. In Fig. || we show the dependence of the Hamming distance D on the system size. As expected from the 
analogy with usual critical phenomena the spreading transition becomes sharper with increasing system size. Figure 
H also shows that a system size of L = 20 already gives reasonable results for the determination of the spreading 
temperature of the Coulomb glass model, provided the disorder strength is comparatively small. 



V. CONCLUSIONS AND OUTLOOK 



We have used the damage-spreading technique to examine the low-temperature dynamics of disordered electronic 
systems with localized states based on the Coulomb glass model. We have found that the dynamics of the system freezes 
below a spreading temperature T s . For zero disorder this damage spreading transition coincides with the equilibrium 
phase transition within our accuracy. At finite disorder strength, when there is no equilibrium phase transition, the 
spreading point T s is shifted to lower temperatures. However, T s remains finite even for larger disorder strengths. 
Consequently, there is a low temperature "phase" of the Coulomb glass with frozen dynamics and a high temperature 
phase where the damage spreads through the system. In the case of finite disorder Wq the spreading transition is 
a purely dynamic transition which does not possess an equilibrium counterpart. A more detailed investigation of 
this transition is in progress. It is, however, hampered by finite-size effects since the long-range interaction severely 
restricts the possible system sizes in our simulations. These limited system sizes are also the reason why the spreading 
point T s for high values of disorder could not yet be determined exactly. 

For small disorder strengths the spreading point T s is still close to the (second-order) equilibrium phase transition 
temperature T c of the system without disorder. Since physical quantities in the vicinity of a critical point can usually 
be described by scaling laws we expect the Hamming distance D to obey the homogeneity relation 

D{W ,T)=tPf(jf), t=\T-T c \ (8) 

with the critical exponents tp and f3. The confirmation of this scaling law and the determination of the exponents 
remain a task for the future. 

One might also ask, how the results change if more sophisticated dynamical algorithms are used, that represent the 
physical processes in disordered insulators better than the simple Metropolis algorithm with single-particle exchange 
with a reservoir. The question is of particular importance, since the properties of damage spreading depend on 
the type of dynamics used in the simulation more strongly than the thermodynamic quantities. We have begun to 
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study the Coulomb glass model with distance-dependent tunneling probabilities between the sites. Results of this 
numerically much more involved investigation will be published elsewhere. 

This work was supported in part by the DA AD, by the DFG under grant number Vo 659/1-1 and SFB 393 and by 
the NSF under grant number DMR-95-10185. 
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